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We report the first observation of seasonal modulations in the rates of cosmic ray multiple-muon 
events at two underground sites, the MINOS Near Detector with an overburden of 225 mwe, 
and the MINOS Far Detector site at 2100 mwe. At the deeper site, multiple-muon events with 
muons separated by more than 8 m exhibit a seasonal rate that peaks during the summer, 
similar to that of single-muon events. In contrast and unexpectedly, the rate of multiple- 
muon events with muons separated by less than 5-8 m, and the rate of multiple-muon events 
in the smaller, shallower Near Detector, exhibit a seasonal rate modulation that peaks in the winter. 
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I. INTRODUCTION 


Muons observed in underground particle detectors 
originate from the interactions of cosmic rays with nu¬ 
clei in the upper atmosphere. These interactions produce 
pions (jt) and kaons ( K ) which can either interact, gen¬ 
erating hadronic cascades, or decay, producing muons. 
The probability that these mesons will decay rather than 
interact is dependent on their energy and the density 
of the atmosphere near their point of production. The 
temperature of the upper atmosphere varies slowly over 
the year, causing a seasonal effect on underground muon 
rates. Increases in the temperature of the atmosphere 
decrease the local density and thus reduce the probabil¬ 
ity that a secondary meson will interact. Consequently, 
the muon flux should increase in the summer. A number 
of experiments have observed this variation in the single 
muon rate [1-11], including MINOS in both Far Detector 
(FD) data fuJliS and Near Detector (ND) data [bj ]. 

Seasonal variations for single muons have been studied 
with a correlation coefficient defined by: 


A Ry 

<R^l> 


A T ef f 

< T eff 


(1) 


such that cxt,n = N x ar,i- The result presented here 
differs greatly from this. This paper presents the first 
measurement of the multiple-muon modulation parame¬ 
ters. 

Note that most extensive air showers have many muons 
in them, but that the highest energy muons which can 
reach an underground detector are produced in the first 
few interactions. Observed single-muon events are most 
likely multiple muons in which any other muons range out 
before reaching the detector or missed the detector later¬ 
ally. A single muon observed in a detector underground 
is most likely the highest energy muon from the shower 
due to the steeply falling cosmic ray energy spectrum. 

The MINOS detectors and the event selection are de¬ 
scribed in Sec. eh In Sec. m the measurement and com¬ 
parison of the modulation parameters for the MINOS ND 
and FD multiple-muon and single-muon event rates are 
presented. In Sec. m and Sec. 0 some possible explana¬ 
tions of the seasonal behavior of the multiple-muon rates 
are considered. 


II. THE MINOS DETECTORS AND MUON 
DATA 


where < R M > is the mean muon rate, and is equivalent to 
the rate for an effective atmospheric temperature equal 
to <T e ff>. The magnitude of the temperature coeffi¬ 
cient (Xt is dependent on the muon energy at production 
and hence the depth of the detector. The effective tem¬ 
perature T e ff is a weighted average over the region of the 
atmosphere where the muons originate. 

By the same reasoning as above a variation should also 
be present in the rate of multiple-muon events. No such 
studies of multiple-muon seasonal rates are reported in 
the literature. The formulae used to calculate T e ff for 
single muons assume a single leading hadron from the 
first interaction is the parent, an assumption that is not 
applicable for multiple-muon events. 

The probability that a cosmic ray shower will give a 
multiple-muon event observed in the MINOS Near or Far 
detectors is enhanced whenever any of the following con¬ 
ditions are true: 1) The primary interaction occurs high 
in the atmosphere where the density is lower and a larger 
fraction of produced hadrons decay; 2) The energy of 
the primary is large so a higher multiplicity of hadrons 
is produced; 3) The cosmic ray primary is a heavy nu¬ 
cleus which breaks up and makes more hadrons; and 4) 
A leading hadron decays to dimuons. Assuming the rel¬ 
ative probability of interaction and decay for each meson 
in a shower is independent, for multiple muons that come 
from the same energy and altitude as a single muon event, 
one might expect an increase in rate during the summer 
that is roughly proportional to the muon multiplicity, N, 


*Now at South Dakota School of Mines and Technology, Rapid 
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The MINOS detectors are planar magnetized 
steel/scintillator tracking calorimeters [l5|. The 
vertically oriented detector planes are composed of 
2.54 cm thick steel and 1cm thick plastic scintillator. 
A scintillator layer is composed of 4.1cm wide strips. 
The MINOS ND has a total mass of 0.98 kton, and lies 
104 m (225 mwe) underground at Fermilab at 42° North 
latitude. The detector is made from 3.8 m x 4.8 m 
hexagonal planes and is 17 m long. It consists of two 
sections, a calorimeter encompassing the upstream 121 
planes and a spectrometer containing the downstream 
161 planes. In both sections, one out of every five planes 
is covered with 96 scintillator strips attached to the steel 
planes. In the calorimeter section, the other four out of 
five planes are covered with 64 scintillator strips, while 
in the spectrometer section they have no scintillator. 
Only muons which enter the calorimeter are included in 
this analysis. The larger FD is 705 m below the surface 
(2100 mwe), has a total mass of 5.4kton, and is located 
in the Soudan Underground Laboratory, at 48° North 
latitude. It is composed of 484 steel-scintillator 8.0 m 
octagonal planes and is 31 m long. The detectors are 
oriented to face the NuMI beam, but through-going 
cosmic muons are well reconstructed over wide geometric 
angular regions. 

Six years of MINOS ND data collected between 
June 1, 2006 and April 30, 2012 and 9years of MINOS FD 
data collected between August 1, 2003 and April 30, 2012 
are analyzed for this paper. The cosmic muon trigger cri¬ 
teria are similar at both detectors requiring that a signal 
is registered in either 4 strips in 5 sequential planes or 
that strips from any 20 planes register a total signal above 
threshold within a given time window. The raw cosmic 
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trigger rate at the ND and FD are approximately 27 Hz 
and 0.5 Hz respectively. 

The single-muon event selection requires there to be a 
single reconstructed track in an event. The multiple- 
muon event selection requires there to be more than 
one reconstructed track in an event. However, since the 
single-muon event rate is much larger than the multiple- 
muon event rate, the multiple-muon sample contains a 
background of single-muon events that have been mis- 
reconstructed to contain two tracks This background is 
greatly reduced by requiring that, for multi-track events, 
the track separation, A S, defined as the minimum point 
of closest approach between any two tracks, be greater 
than 0.6 m. Observed excesses due to this background 
at small AS in both the ND and FD were removed by 
this selection, reducing the background from 1.3% to less 
than 200 events out of 11 million in the FD. Figure [1] 
shows the time between sequential multiple-muon events. 
The multiple-muon event rates at the ND and FD are 
19.6 mHz and 14.1 mHz respectively. In total the MI¬ 
NOS ND and FD have collected 2.45 xlO 6 and 3.36 xlO 6 
good multiple-muon events respectively. 

The rate of multiple muons in the MINOS detectors is 
dominated by m^ = 2 and m M = 3, where is the muon 
multiplicity. For the FD, the reconstruction works well 
for these small multiplicities, identifying all tracks in 84% 
(67%) of events for m M = 2 (3). From a multiple-muon 
Monte Carlo [il|, the efficiency for identifying an event 
as a multiple muon is 84% (94%) for m M = 2 (3), rising 
to above 97% for m M > 3. However, the reconstructed 
multiplicity is frequently too low for high multiplicity 
events. No event with > 13 is recorded in either 
of the MINOS detectors. Similar measurements with the 
finer-grained Soudan 2 detector at the same depth as the 
MINOS FD recorded multiplicities up to 20 (n|- In the 
coarser-grained MINOS detectors, the individual tracks 
closest in distance from such high-multiplicity events will 
be resolved as a single track or not pass the track quality 
criteria used in the MINOS reconstruction algorithms. 
Figure [2] shows the reconstructed multiplicity distribu¬ 
tion in the MINOS FD. 


III. MODULATION ANALYSIS 

To compare the variation in the event rates for 
multiple-muon and single-muon events, the rates are fit 
to a sinusoidally-varying function of time. There is no 
a-priori reason to believe that the rates vary sinusoidally 
through the year, but this fit gives a qualitatively use¬ 
ful amplitude and phase. The following function, which 
contains four free parameters, is used for the fit: 

m = Ml ~ 3 ^g)(l + Acos[^(t - to)]) (2) 

where t is the number of days since Jan. 1, 2010 and f 0 
is the phase; Rg is the mean rate on Jan. 1, 2010; A 



FIG. 1: Time between neighboring atmospheric multiple- 
muon events in the MINOS detectors. The data are well 
described by an exponential over six orders of magnitude in 
instantaneous rate. 



FIG. 2: The reconstructed muon multiplicity, for events con¬ 
taining more than one reconstructed track, in the Far Detec¬ 
tor. 


is the modulation amplitude and T is the period (ap¬ 
proximately 1 year). The parameter / is the loss rate 
(described in Reference [lj]) that accounts for an ob¬ 
served linear decrease in the event rate in both the FD 
and ND over the lifetime of the experiment. The source 
of this small but apparently steady decrease has not been 
conclusively identified and is under study. The best-fit 
parameters are given in Table |T] 
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Data Set 

Amplitude 

Loss Rate (f) Period (T) 

Phase (to) 


(%) 

(%/year) 

(days) 

(days) 

MINOS FD 

AS>0.6m 

0.39T0.08 

-0.04i0.02 

356.4T4.1 

105.2T16.1 

0.6 m < AS < 4.5 m 

l.OiO.l 

-0.14i0.04 

362.2F3.3 

27.6i8.9 

4.5 m < AS < 8.0 m 

0.47T0.14 

0.02i0.04 

354.6i9.1 

78.9il7.3 

AS > 8.0m 

2.0T0.1 

0.01T0.04 

363.7il.8 

184.8T6.5 

Single Muons 

1.27T0.01 

0.013T0.001 

364.4T0.3 

183.0T0.9 

MINOS ND 

AS > 0.6 m 

2.51T0.09 

0.35T0.03 

367.4il.3 

23.7A2.3 

0.6 m < AS < 1.8 m 

2.35T0.17 

0.25T0.05 

369.0A2.5 

26.2A4.2 

1.8 m < AS < 3.0 m 

2.53T0.17 

0.41T0.05 

369.3A2.3 

25.1T4.0 

AS > 3.0 m 

2.64T0.17 

0.39T0.05 

365.8A2.1 

22.1A3.8 

Single Muons 

0.268A0.004 

0.0116A0.001 

365.7T0.4 

198.6T0.9 


TABLE I: The parameters obtained when Eq. [2] is fit to the single-muon and multiple-muon data in each detector. The table 
also shows the results of fits to subsets of the multiple-muon data, based on the minimum separation between tracks. The best 
fit phase and period do not change significantly if the loss rate is assumed to be zero. 


A. Modulations in the Far Detector 

The fit for seasonal variations in the FD multiple-muon 
sample shows a much smaller amplitude than for single 
muons, and a poorly defined phase. Since the MINOS 
FD is larger than the ND and is fully instrumented, the 
modulation is studied as a function of track separation. 
Figure [3] shows the track separation AS. The multiple- 
muon data are grouped into three bins of roughly equal 
statistics with track separations from 0.6-4.5 m (FD re¬ 
gion A), 4.5-8.0m (FD region B) and greater than 8m 
(FD region C). Region A most closely resembles the dis¬ 
tribution in the ND. 

Figure 0] presents the multiple-muon rate in the MI¬ 
NOS FD as a function of time for differing track separa¬ 
tions. The FD multiple-muon data set with the largest 
track separation, >8 m, modulates with a summer maxi¬ 
mum (to = 184.8±6.5days); this phase is consistent with 
that observed in the FD single-muon sample, and the am¬ 
plitude is larger. On the other hand, the FD multiple- 
muon data set with the smallest track separations mod¬ 
ulates with a winter maximum (to = 27.6±8.9days); 

this phase differs by a half year from the variation seen 
with single muons. The FD mid-range track separation 
multiple-muon data set has a small amplitude and is con¬ 
sistent with an admixture of the other two phases. 

In Fig. 0 the data for regions A and C have been 
binned by calendar month, with each point showing the 
average rate over all years of data-taking. 


B. Modulations in the Near Detector 

The ND multiple-muon data, shown in Fig. [Gl and the 
single-muon data (shown in Reference 0 ) were fit to 
Eq. fusing one month time interval bins. The multiple- 



FIG. 3: The minimum track separation AS between any 
two tracks in multiple-muon events recorded in the FD. The 
gray (black) histogram is the distribution before (after) the 
selection to remove misreconstructed single-muon events. Re¬ 
gions of track separation AS are defined as A: 0.6-4.5 m, B: 
4.5-8.0m and C: > 8 m. 


muon event rate data show a clear modulation signature. 
However, unlike the single muon rate which reaches its 
maximum in the summer 0, the multiple-muon rate 
reaches its maximum in the winter. This also matches 
the modulation for the region-A multiple muons in the 
FD. Both the single-muon and multiple-muon data sets 
have periods consistent with one year but their phases, 
198.6±0.9days and 23.7±2.3days respectively, differ by 
about six months. The rates of multiple muons and single 
muons, binned by calendar month and averaged over all 
years of data-taking, are shown in Fig. [7] 
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FIG. 4: The multiple-muon rate in the FD as a function of 
time for different track separations. Each data point corre¬ 
sponds to one calendar month of data. The solid red lines are 
the best fit to Eq. [2] The top graph is for the smallest track 
separation, the middle graph for mid-range and the bottom 
graph for the largest. The vertical lines are year boundaries 
and the solid horizontal line represents the fit without the 
cosine term. 

Figure [ 8 ] shows the track separation in ND multiple- 
muon events. To qualitatively match the procedure in the 
FD, the data have been grouped into three bins of roughly 
equal statistics with track separations of 0.6-1 .8 m (ND 
region A), 1.8-3.0m (ND region B) and greater than 3m 
(ND region C). As before, the data are fit to Eq. d2]) 
and the best fit parameters are given in Table U| There 
is no apparent difference in the fit parameters for the 
three ND regions, which all peak in the winter. There is 
consistency between ND regions ABC and FD region A 
in both AS and a winter maximum. 


IV. DISCUSSION OF RESULTS AND POSSIBLE 
EXPLANATIONS 

We have previously observed seasonal variations in 
single-muon rates in the MINOS ND and FD that corre¬ 
late at expected levels with the temperature changes and 
the season. Those muon rates rose in the summer as did 
the calculated values of T e ff , and the measured correla¬ 
tions were q™ = 0.428 ± 0.059 [HJ and q™ = 0.873 
± 0.014 [l3j]- The measurement of a multiple-muon rate 
in the ND that peaks in the winter is unexpected, as is 
the winter maximum in the FD in region A of separation. 
In order to try to understand this result, four plausible 
explanations which might account for these results are 



FIG. 5: The multiple-muon rate in the FD for events with 
AS range A from 0.6 m to 4.5 m (top graph) and for events 
with AS range C larger than 8 m (bottom) binned according 
to calendar month. The top figure shows a winter maximum. 
The bottom figure shows a summer maximum. 
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FIG. 6: The multiple-muon rate in the ND as a function of 
time. Each data point corresponds to one calendar month. A 
clear modulation in the data is observed with the maximum 
occurring towards the start of the year. The vertical lines are 
year boundaries. 


considered. They involve: A) a source of dimuons from 
prompt hadron decays (such as 77 and p) that may have 
the opposite seasonal variation, since in the winter the 
secondary pions are more likely to interact than decay 
and produce more of such hadrons; B) a geometric effect 
in which different altitude distributions affect the track 
separation underground; C) a different altitude distribu¬ 
tion for multimuon events that may come from regions 
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FIG. 8: The minimum track separation AS between any 
two tracks in multiple-muon events recorded in the ND. The 
gray (black) histogram is the distribution before(after) the se¬ 
lection to remove misreconstructed single-muon events. Re¬ 
gions of track separation AS are defined as A: 0.6-1.8 m, B: 
1.8-3.0m and C: > 3 m. 


ability of other hadrons which have dimuon decays, such 
as ?? and p mesons, changes negligibly. The 2% more pi- 
ons [3 which interact will increase the number of these 
other hadrons. This increase, which is at most 2%, must 
then be folded in with the small dimuon branching ra¬ 
tios, such as 4.6 x 10 -5 for p —>• p + p~ and 3.1 x 10 “ 4 
for 77 —> p + p ~7 [35j. Observed dimuon rates are 1% of 
the single muon rates in the FD, and 0.16% in the ND, 
so even if p and 77 production were comparable to 7r, this 
contribution is at most 6 x 10 -6 , too small to account 
for the observed effect. 


FIG. 7: The top figure is the multiple-muon rate in the 
ND, binned according to calendar month, which each point 
showing the average rate for all years of data-taking. The 
figure also shows a cosine fit to the data. The single-muon 
rate is shown in the bottom figure, showing a clearly different 
seasonal modulation. 


of the atmosphere with different seasonal temperature 
profiles; and D) leading secondary hadrons being more 
likely to decay than interact in the summer, and thus 
less likely to make multiple hadrons which make multi¬ 
ple muons. We discuss each of these possibilities in the 
current section. 


A. Hadronic dimuon decays 

One idea is that the winter maximum may be due to 
hadronic decays into dimuons. In the winter, while pions 
are less likely to decay in the atmosphere, the decay prob- 


B. A geometry effect 

A possibility is that the muons generated higher in the 
atmosphere in the summer spread out farther so that 
there are fewer of them in region A. This would be solely 
a geometric effect, in that it would not affect the num¬ 
ber of multimuons in each season but only the track- 
separation distribution. This is further complicated by 
multiple scattering, but an effect due to the opening an¬ 
gle at production can be estimated. For a fixed-size de¬ 
tector, a difference in the track separation distribution 
would affect the measured rate. The altitude of the first 
interaction in an isothermal atmosphere is related to the 
absolute temperature. A ± 2% seasonal change in the 
effective temperature would cause a ± 2 % change in the 
altitude, and hence less than a 4% change in the average 
muon track separation underground. This would move 
events to the right in Fig. [3] Due to the shape of the 
distribution, more events would move from region A to 
region B than from region B to region C, which is in 
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variation at those altitudes. 

In order to study the possible altitude dependence 
of multiple muons we simulated cosmic ray air showers 
which could make multiple muons in the MINOS FD. The 
Monte Carlo sample was produced by CORSIKA uni 
using version 7.4. We have run CORSIKA with three 
different hadronic models, QGSJET-01C, QGSJET 11-04 
[32| and EPOS [33| which gave consistent results. We 
note that CORSIKA uses an isothermal atmosphere and 
cannot be used per se to study seasonal variations. The 
goal here is to roughly calculate the altitude dependence 
for the three regions of track separation. CORSIKA out¬ 
puts muon energies and positions at the earth’s surface. 
To reach the MINOS FD, energy loss through the rock 
was calculated using [34j |: 



CD 
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FIG. 9: The (top) modulation phase and (bottom) ampli¬ 
tude in the ECMWF temperature data based on a cosine fit 
are shown as a function of altitude and detector site. These 
distributions were used to study both the geometry effect (B) 
and the temperature effect (C). 


contradiction to our fits. Also, one would expect a sim¬ 
ilar effect in the ND as shown in Fig. [Si but no track 
separation dependence is seen in the ND. 


C. A temperature effect 


Eioss(X) = ^-{e bTX - 1), (3) 

where X is the rock overburden, a is a parameter for the 
ionization energy loss and = Hrem + b pa i r + bois rep¬ 
resents the energy loss due to bremsstrahlung, electron- 
positron pair production and photo-nuclear interactions. 

Simulated events were selected for which two or more 
muons reached the top of the FD with a total remaining 
energy of at least 0.9 GeV. The distribution of track sep¬ 
aration obtained with this simulation was similar to, but 
not identical to, the distribution seen in data (Fig. [3]). 
We then extracted from CORSIKA the altitude at which 
each muon was created in each track-separation region. 
Those three distributions are shown in Fig.QjI] There is a 
shift in the mean altitude for each region of track separa¬ 
tion from 17 km in region A to 21 km in region C, though 
all three distributions are quite broad. We then combined 
the altitude dependence with the temperature phase and 
amplitude fits shown in Fig. [9l assuming the rate and 
T e ff were completely correlated, to compare the over¬ 
all variation of T e f f averaged over each track-separation 
region. The result was a variation that peaked in the 
summer in all three regions, with an amplitude of 1.9% in 
region C and 1.6% in regions A and B. This study was re¬ 
peated using QGSJET-01C, QGSJET II and EPOS and 
all three results were similar. It does not appear that the 
temperature variations noted in Fig. [9] can account for 
the observed reverse seasonal effect in region A. 


To determine whether there may be an altitude- 
dependent seasonal variation that differs for single and 
multiple muons, meteorological data is used to deter¬ 
mine the atmospheric temperature profile. Figure [9] gives 
the phase and amplitude of the modulation of the atmo¬ 
spheric temperature, based on a cosine fit to data taken 
from the European Center for Medium-Range Weather 
Forecasts (ECMWF) model [29|, as a function of at¬ 
mospheric pressure. Indeed, there is a small region of 
the atmosphere, between 70hPa and 175 hPa, where the 
temperature reaches a maximum in the winter. Note, 
however, the small amplitude of the annual temperature 


D. Anticorrelation of primary and secondary 
decays 

As a last hypothesis, while most single-muon events 
come from secondary pions and kaons produced in the 
primary cosmic ray interaction, multiple muons may be 
more likely to come from higher energy primaries where 
there are further hadronic interactions deeper in the 
shower. In that case, if the secondary hadron is more 
likely to decay in the summer, it is less likely to interact 
and make additional pions and kaons which contribute to 
multiple muons. This may be the best explanation for the 





FIG. 10: To study a possible temperature effect with altitude, 
fSec. IIV~Cl in the text), the altitude distribution from COR- 
SIKA for MINOS FD multiple muons are shown for each of 
the three regions of track separation in Fig. [3] 

winter maximum measured in the MINOS ND multiple- 
muon data set. A quantitative test of this hypothesis 
will require a detailed study of air shower development 
that is beyond the scope of this analysis. This hypoth¬ 
esis accounts for the stronger effect in the MINOS ND, 
where the muons come from pions and kaons below their 
critical energies (e n = 115 GeV and ex = 850 GeV, de¬ 
fined as those energies for which meson decay and inter¬ 
action rates in the atmosphere where muons originate are 
equal) [l] and for the more complex effect in the MINOS 
FD where the energies are above and comparable to 
ex ■ Mesons which are much below their critical energies 
mostly decay, so the temperature effect that does exist to 
increase the decay rate in the summer has a large effect 
on decreasing the interaction rate in the summer. This 
is the situation for muons in the ND where the threshold 
from the overburden is near 50 GeV. At the FD, where 
the threshold is almost a TeV, a change in the decay rate 
has a smaller impact on the interaction rate, since a large 
fraction of the hadrons are interacting before they decay. 

As pointed out in the introduction, single muons come 
predominantly from the decay of a leading hadron, and 


multiple muons from a more complicated process. It is 
clear that if a leading hadron is more likely to decay in 
one season, it is less likely to interact. 


V. CONCLUSION 

We have shown evidence of an annual modulation in 
the MINOS ND multiple-muon data set in which the 
maximum rate occurs in the winter. This phase is incon¬ 
sistent with the summer maximum observed in the ND 
and FD single-muon data. Data collected by the MINOS 
FD were used to show that there is a transition from a 
summer maximum in multiple-muon events with a large 
track separation to a winter maximum in multiple-muon 
events with a small track separation. This transition oc¬ 
curs at track separations of about 5-8 m. 

Four possible explanations for this observed character¬ 
istic in seasonal variations were considered. One expla¬ 
nation is favored: this is a hypothesis in which multiple 
muons come preferentially from higher energy pions and 
kaons which, in the summer, are less likely to interact and 
produce the secondary pions and kaons that give rise to 
the multiple muons. However, a full explanation of our 
observations including the dependence in the FD on track 
separation must come from a more detailed study of ex¬ 
tensive air-shower properties and the properties of the 
atmosphere. 
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